R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

library(raster)
size_dat=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/fire_growth_5days_v4.txt", header=T, row.names=NULL)

size_dat=as.data.frame(size_dat)
size_dat=size_dat[-1,]
size_dat=size_dat[,-1]
colnames(size_dat)=c("firename","year","cause","size1","size2","size3","size4","size5","final_firesize","mean_precip1","mean_precip2","mean_precip3","mean_precip4","mean_precip5","mean_tmax1","mean_tmax2","mean_tmax3","mean_tmax4","mean_tmax5","mean_tmean1","mean_tmean2","mean_tmean3","mean_tmean4","mean_tmean5","mean_vpdmax1","mean_vpdmax2","mean_vpdmax3","mean_vpdmax4","mean_vpdmax5","mean_windspeed1","mean_windspeed2","mean_windspeed3","mean_windspeed4","mean_windspeed5","landcover","ecosystem","biomass","elevation")

size_dat2 <- data.frame(lapply(size_dat, function(x) as.numeric(as.character(x))))

NAs introduced by coercion

size_dat2$human = 0
size_dat2$human[size_dat2$cause !=1 & size_dat2$cause !=14 & size_dat2$cause !=17]=1
size_dat2$human[size_dat2$cause ==1 ]=2

pro1 =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 1 ), ]
pro2 =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 1 ), ]

t.test(pro1$size1,pro2$size1)

    Welch Two Sample t-test

data:  pro1$size1 and pro2$size1
t = -2.0397, df = 74.381, p-value = 0.04493
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -22.0986230  -0.2593373
sample estimates:
mean of x mean of y 
 3.056385 14.235365 

t.test(pro1$size2,pro2$size2)

    Welch Two Sample t-test

data:  pro1$size2 and pro2$size2
t = -2.758, df = 73.138, p-value = 0.007341
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -32.702045  -5.266208
sample estimates:
mean of x mean of y 
 8.254394 27.238521 

t.test(pro1$size3,pro2$size3)

    Welch Two Sample t-test

data:  pro1$size3 and pro2$size3
t = -2.9002, df = 58.203, p-value = 0.005254
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -51.965993  -9.526741
sample estimates:
mean of x mean of y 
 14.10484  44.85121 

t.test(pro1$size4,pro2$size4)

    Welch Two Sample t-test

data:  pro1$size4 and pro2$size4
t = -3.1797, df = 53.078, p-value = 0.002461
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -74.99798 -16.98040
sample estimates:
mean of x mean of y 
 17.68430  63.67349 

t.test(pro1$size5,pro2$size5)

    Welch Two Sample t-test

data:  pro1$size5 and pro2$size5
t = -3.4218, df = 39.443, p-value = 0.001462
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -117.71273  -30.26844
sample estimates:
mean of x mean of y 
 19.71795  93.70853 

pro =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 1), ]
length(pro$year)

[1] 68

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

pro =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 1), ]
length(pro$year)

[1] 77

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

pro =size_dat2[which(size_dat2$human == 1 & size_dat2$landcover == 2), ]
length(pro$year)

[1] 18

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

pro =size_dat2[which(size_dat2$human == 2 & size_dat2$landcover == 2), ]
length(pro$year)

[1] 9

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

length(!is.na(pro$size5))

[1] 9

pro =size_dat2[which(size_dat2$human == 1 & size_dat2$ecosystem == 6), ]
length(pro$year)

[1] 41

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

pro =size_dat2[which(size_dat2$human == 2& size_dat2$ecosystem == 6), ]
length(pro$year)

[1] 79

boxplot(pro$size1,pro$size2,pro$size3,pro$size4,pro$size5,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,500), cex.lab=1.4,cex.axis = 1.3)

NA
NA

plot fire size map QGIS

library(data.table)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
data.table 1.13.0 using 1 threads (see ?getDTthreads).  Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********

Attaching package: ‘data.table’

The following object is masked from ‘package:raster’:

    shift

DT= data.table(res)
fire_size = DT[ , .SD[which.min(growth_km)], by = firename]

Error in .checkTypos(e, names_x) : 
  Object 'growth_km' not found amongst bi, erc, etr, fm100, fm1000 and 23 more

extract number and size statistics from frap

library(raster)

dr1 =shapefile("/Users/stijnhantson/Documents/data/FRAP/fire18_1.shp")
dr1$YEAR_=as.numeric(as.character(dr1$YEAR_))
dr1$Shape_Area=as.numeric(as.character(dr1$Shape_Area))
dr1$ALARM_DATE = as.Date(dr1$ALARM_DATE )
dr1=dr1[!is.na(dr1$YEAR_), ]
dr1=dr1[dr1$YEAR_>2011,]


dr2=dr1[dr1$YEAR_ == 2013,]
plot(table(dr2$ALARM_DATE), xlim=c(as.Date("2013/01/01"),as.Date("2013/12/31")))

table(factor(dr2$ALARM_DATE, levels = as.Date("2013/01/01"):as.Date("2013/12/31")))

15706 15707 15708 15709 15710 15711 15712 15713 15714 15715 15716 15717 15718 15719 15720 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15721 15722 15723 15724 15725 15726 15727 15728 15729 15730 15731 15732 15733 15734 15735 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15736 15737 15738 15739 15740 15741 15742 15743 15744 15745 15746 15747 15748 15749 15750 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15751 15752 15753 15754 15755 15756 15757 15758 15759 15760 15761 15762 15763 15764 15765 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15766 15767 15768 15769 15770 15771 15772 15773 15774 15775 15776 15777 15778 15779 15780 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15781 15782 15783 15784 15785 15786 15787 15788 15789 15790 15791 15792 15793 15794 15795 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15796 15797 15798 15799 15800 15801 15802 15803 15804 15805 15806 15807 15808 15809 15810 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15811 15812 15813 15814 15815 15816 15817 15818 15819 15820 15821 15822 15823 15824 15825 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15826 15827 15828 15829 15830 15831 15832 15833 15834 15835 15836 15837 15838 15839 15840 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15841 15842 15843 15844 15845 15846 15847 15848 15849 15850 15851 15852 15853 15854 15855 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15856 15857 15858 15859 15860 15861 15862 15863 15864 15865 15866 15867 15868 15869 15870 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15871 15872 15873 15874 15875 15876 15877 15878 15879 15880 15881 15882 15883 15884 15885 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15886 15887 15888 15889 15890 15891 15892 15893 15894 15895 15896 15897 15898 15899 15900 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15901 15902 15903 15904 15905 15906 15907 15908 15909 15910 15911 15912 15913 15914 15915 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15916 15917 15918 15919 15920 15921 15922 15923 15924 15925 15926 15927 15928 15929 15930 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15931 15932 15933 15934 15935 15936 15937 15938 15939 15940 15941 15942 15943 15944 15945 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15946 15947 15948 15949 15950 15951 15952 15953 15954 15955 15956 15957 15958 15959 15960 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15961 15962 15963 15964 15965 15966 15967 15968 15969 15970 15971 15972 15973 15974 15975 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15976 15977 15978 15979 15980 15981 15982 15983 15984 15985 15986 15987 15988 15989 15990 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
15991 15992 15993 15994 15995 15996 15997 15998 15999 16000 16001 16002 16003 16004 16005 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
16006 16007 16008 16009 16010 16011 16012 16013 16014 16015 16016 16017 16018 16019 16020 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
16021 16022 16023 16024 16025 16026 16027 16028 16029 16030 16031 16032 16033 16034 16035 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
16036 16037 16038 16039 16040 16041 16042 16043 16044 16045 16046 16047 16048 16049 16050 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
16051 16052 16053 16054 16055 16056 16057 16058 16059 16060 16061 16062 16063 16064 16065 
    0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 
16066 16067 16068 16069 16070 
    0     0     0     0     0 

human lightning unknown ratio through time

library(raster)
library(rgdal)

package ‘rgdal’ was built under R version 3.6.2rgdal: version: 1.5-12, (SVN revision 1018)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
GDAL binary built with GEOS: FALSE 
Loaded PROJ runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
Linking to sp version:1.4-2
Overwritten PROJ_LIB was /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj

library(rgeos)

package ‘rgeos’ was built under R version 3.6.2rgeos version: 0.5-3, (SVN revision 634)
 GEOS runtime version: 3.7.2-CAPI-1.11.2 
 Linking to sp version: 1.4-1 
 Polygon checking: TRUE 

library(sf)
frap=readOGR("/Users/stijnhantson/Documents/data/FRAP/fire19_1.shp")

OGR data source with driver: ESRI Shapefile 
Source: "/Users/stijnhantson/Documents/data/FRAP/fire19_1.shp", layer: "fire19_1"
with 20820 features
It has 18 fields
Integer64 fields read as strings:  CAUSE C_METHOD OBJECTIVE 

shape = shapefile("/Users/stijnhantson/Documents/data/veg_california/ca_eco_l3/ca_eco_l3.shp")
shape = spTransform(shape,crs(frap))


frap = st_make_valid(st_read("/Users/stijnhantson/Documents/data/FRAP/fire19_1.shp"))

Reading layer `fire19_1' from data source `/Users/stijnhantson/Documents/data/FRAP/fire19_1.shp' using driver `ESRI Shapefile'

GDAL Message 1: organizePolygons() received an unexpected geometry.  Either a polygon with interior rings, or a polygon with less than 4 points, or a non-Polygon geometry.  Return arguments as a collection.GDAL Message 1: Geometry of polygon of fid 19691 cannot be translated to Simple Geometry. All polygons will be contained in a multipolygon.

Simple feature collection with 20820 features and 18 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -373237.5 ymin: -604727.6 xmax: 519987.8 ymax: 518283.7
CRS:            3310

types <- vapply(sf::st_geometry(frap), function(x) {
  class(x)[2]
}, "")
polys <- frap[ grepl("*POLYGON", types), ]

frap=sf:::as_Spatial(polys)
frap=gBuffer(frap, byid=T,  width=0.0)
eco_inter=intersect(frap,shape)

NA

1. frap
2. only fire growth datasets

• prepare final dataset to open

library(raster)
#library(rgdal)
daily_res=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/final_dataset_V5.txt",header=T)

res=as.data.frame(daily_res)

res$mean_ros =as.numeric(as.character(res$mean_ros))
res$max_ros =as.numeric(as.character(res$max_ros))
res$median95_ros =as.numeric(as.character(res$median95_ros))
res$bi =as.numeric(as.character(res$bi))
res$erc =as.numeric(as.character(res$erc))
res$etr =as.numeric(as.character(res$etr))
res$fm100 =as.numeric(as.character(res$fm100))
res$fm1000 =as.numeric(as.character(res$fm1000))
res$pet =as.numeric(as.character(res$pet))
res$pr =as.numeric(as.character(res$pr))
res$rmax =as.numeric(as.character(res$rmax))
res$rmin =as.numeric(as.character(res$rmin))
res$th =as.numeric(as.character(res$th))
res$tmmn =as.numeric(as.character(res$tmmn))
res$tmmx =as.numeric(as.character(res$tmmx))
res$vpd =as.numeric(as.character(res$vpd))
#res$ws =as.numeric(as.character(res$ws))
res$vs =as.numeric(as.character(res$vs))
res$growth =as.numeric(as.character(res$growth))
res$total_area =as.numeric(as.character(res$total_area))
res$mean_frp =as.numeric(as.character(res$mean_frp))
res$frp_95 =as.numeric(as.character(res$frp_95))
res$max_land =as.numeric(as.character(res$max_land))
res$mean_land =as.numeric(as.character(res$mean_land))
res$biomass =as.numeric(as.character(res$biomass))
res$year =as.numeric(as.character(res$year))
res$month =as.numeric(as.character(res$month))
res$doy_out =as.numeric(as.character(res$doy_out))

res$mean_dnbr =as.numeric(as.character(res$mean_dnbr))
res$mean_rdnbr =as.numeric(as.character(res$mean_rdnbr))
res$mean_bas =as.numeric(as.character(res$mean_bas))
res$median_dnbr =as.numeric(as.character(res$median_dnbr))
res$median_rdnbr =as.numeric(as.character(res$median_rdnbr))
res$median_bas =as.numeric(as.character(res$median_bas))
res$q95_dnbr =as.numeric(as.character(res$q95_dnbr))
res$q95_rdnbr =as.numeric(as.character(res$q95_rdnbr))
res$q95_bas =as.numeric(as.character(res$q95_bas))


res = res[-1,]
res$per_ba = res$growth/res$total_area
res$growth_km =res$growth/1000000

res$human = 0
res$human[res$cause !=1 & res$cause !=14 & res$cause !=17]=1
res$human[res$cause ==1 ]=2

res$ros_km = (res$median95_ros*24)/1000
res$ros_mean_km = (res$mean_ros*24)/1000

is there relation between days untill 75% and fire size

dr1 =shapefile("/Users/stijnhantson/Documents/data/FRAP/fire18_1.shp")
dr1$YEAR_=as.numeric(as.character(dr1$YEAR_))
dr1$Shape_Area=as.numeric(as.character(dr1$Shape_Area))
dr1=dr1[!is.na(dr1$YEAR_), ]
dr1=dr1[dr1$YEAR_>2011,]

fi = dr1[dr1$OBJECTID==17495,]
dr1=subset(dr1, GIS_ACRES>300)

summary(dr1)

Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x -332881.6 349023.3
y -601475.0 462306.3
Is projected: TRUE 
proj4string :
[+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs]
Data attributes:
    OBJECTID         YEAR_         STATE              AGENCY            UNIT_ID           FIRE_NAME           INC_NUM         
 Min.   :11970   Min.   :2012   Length:515         Length:515         Length:515         Length:515         Length:515        
 1st Qu.:17777   1st Qu.:2014   Class :character   Class :character   Class :character   Class :character   Class :character  
 Median :19129   Median :2016   Mode  :character   Mode  :character   Mode  :character   Mode  :character   Mode  :character  
 Mean   :19125   Mean   :2015                                                                                                 
 3rd Qu.:20512   3rd Qu.:2017                                                                                                 
 Max.   :21127   Max.   :2018                                                                                                 
                                                                                                                              
  ALARM_DATE         CONT_DATE            CAUSE             COMMENTS           REPORT_AC          GIS_ACRES       
 Length:515         Length:515         Length:515         Length:515         Min.   :     0.0   Min.   :   303.8  
 Class :character   Class :character   Class :character   Class :character   1st Qu.:   504.1   1st Qu.:   612.4  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character   Median :  1651.5   Median :  1707.6  
                                                                             Mean   : 12021.3   Mean   : 12138.1  
                                                                             3rd Qu.:  6970.0   3rd Qu.:  7068.0  
                                                                             Max.   :410203.0   Max.   :410202.5  
                                                                             NA's   :67                           
   C_METHOD          OBJECTIVE           FIRE_NUM           Shape_Leng       Shape_Area       
 Length:515         Length:515         Length:515         Min.   :  4635   Min.   :1.229e+06  
 Class :character   Class :character   Class :character   1st Qu.: 10017   1st Qu.:2.478e+06  
 Mode  :character   Mode  :character   Mode  :character   Median : 17475   Median :6.910e+06  
                                                          Mean   : 41493   Mean   :4.912e+07  
                                                          3rd Qu.: 37905   3rd Qu.:2.860e+07  
                                                          Max.   :445282   Max.   :1.660e+09  
                                                                                              

res75 = res[res$per_ba > 0.999,]

peak_day1 = as.data.frame(aggregate(res75$fire_day, by = list(res75$firename,res75$year), min))
p=13
year1=0
first=0
last=0
len_d=0
day75=0
firenam1=0
name_fire=0
size=0
k=0
for (p in 1:(length(peak_day1$Group.1)))
{
  print(p)
  year = peak_day1[p,2]
    firenam = as.character(peak_day1[p,1]) 
  
  fir = dr1[which(dr1$YEAR_==year & dr1$FIRE_NAME==firenam),]
  if (length(fir)>1){
    maxsize = max(fir$GIS_ACRES)
    fir = fir[fir$GIS_ACRES ==  maxsize,]
  }
  
  if (!is.na(fir$ALARM_DATE)  | !is.na(fir$ALARM_DATE)){
    k=k+1
    day75[k] = peak_day1[p,3]
    firenam1[k]=firenam 
  year1[k] = year
  first[k]= fir$ALARM_DATE
  last[k] = fir$CONT_DATE
  len_d[k] = as.numeric(as.Date(fir$CONT_DATE) -as.Date(fir$ALARM_DATE))
 name_fire[k]= fir$OBJECTID
    size[k] = fir$GIS_ACRES
  }
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50
[1] 51
[1] 52
[1] 53
[1] 54
[1] 55
[1] 56
[1] 57
[1] 58
[1] 59
[1] 60
[1] 61
[1] 62
[1] 63
[1] 64
[1] 65
[1] 66
[1] 67
[1] 68
[1] 69
[1] 70
[1] 71
[1] 72
[1] 73
[1] 74
[1] 75
[1] 76
[1] 77
[1] 78
[1] 79
[1] 80
[1] 81
[1] 82
[1] 83
[1] 84
[1] 85
[1] 86
[1] 87
[1] 88
[1] 89
[1] 90
[1] 91
[1] 92
[1] 93
[1] 94
[1] 95
[1] 96
[1] 97
[1] 98
[1] 99
[1] 100
[1] 101
[1] 102
[1] 103
[1] 104
[1] 105
[1] 106
[1] 107
[1] 108
[1] 109
[1] 110
[1] 111
[1] 112
[1] 113
[1] 114
[1] 115
[1] 116
[1] 117
[1] 118
[1] 119
[1] 120
[1] 121
[1] 122
[1] 123
[1] 124
[1] 125
[1] 126
[1] 127
[1] 128
[1] 129
[1] 130
[1] 131
[1] 132
[1] 133
[1] 134
[1] 135
[1] 136
[1] 137
[1] 138
[1] 139
[1] 140
[1] 141
[1] 142
[1] 143
[1] 144
[1] 145
[1] 146
[1] 147
[1] 148
[1] 149
[1] 150
[1] 151
[1] 152
[1] 153
[1] 154
[1] 155
[1] 156
[1] 157
[1] 158
[1] 159
[1] 160
[1] 161
[1] 162
[1] 163
[1] 164
[1] 165
[1] 166
[1] 167
[1] 168
[1] 169
[1] 170
[1] 171
[1] 172
[1] 173
[1] 174
[1] 175
[1] 176
[1] 177
[1] 178
[1] 179
[1] 180
[1] 181
[1] 182
[1] 183
[1] 184
[1] 185
[1] 186
[1] 187
[1] 188
[1] 189
[1] 190
[1] 191
[1] 192
[1] 193
[1] 194
[1] 195
[1] 196
[1] 197
[1] 198
[1] 199
[1] 200
[1] 201
[1] 202
[1] 203
[1] 204
[1] 205
[1] 206
[1] 207
[1] 208
[1] 209
[1] 210
[1] 211
[1] 212

plot(len_d,day75, xlab="fire duration FRAP (days)", ylab="fire duration VIIRS (days)")

plot mean vs max fire rate-of-spread

#summary(res)
plot(res$ros_km,res$ros_mean_km, xlab="maximum fire-spread-rate (km/day",ylab="mean fire-spread-rate (km/day)")


pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/mean_vs_max_ros_v1.pdf", width = 7, height = 7)

plot(res$ros_km,res$ros_mean_km, xlim=c(0,25),ylim=c(0,10), xlab="fire rate-of-spread (km/day)",ylab="mean fire rate-of-spread (km/day)", cex.lab=1.3,cex.axis = 1.25)
dev.off()

quartz_off_screen 
                2 

difference between human and lightnign fires

me=0
me1=0
days = c("day1","day2","day3","day4","day5")
days1 = c("1","2","3","4","5")
pro1 = res[res$fire_day == 1 & res$human == 1,]
pro2 = res[res$fire_day == 2 & res$human == 1,]
pro3 = res[res$fire_day == 3 & res$human == 1,]
pro4 = res[res$fire_day == 4 & res$human == 1,]
pro5 = res[res$fire_day == 5 & res$human == 1,]
me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

pro1h = res[res$fire_day == 1 & res$human == 2,]
pro2h = res[res$fire_day == 2 & res$human == 2,]
pro3h = res[res$fire_day == 3 & res$human == 2,]
pro4h = res[res$fire_day == 4 & res$human == 2,]
pro5h = res[res$fire_day == 5 & res$human == 2,]
me1[1] =mean(pro1h$growth_km,na.omit=T)
me1[2] =mean(pro2h$growth_km,na.omit=T)
me1[3] =mean(pro3h$growth_km,na.omit=T)
me1[4] =mean(pro4h$growth_km,na.omit=T)
me1[5] =mean(pro5h$growth_km,na.omit=T)

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,250), cex.lab=1.4,cex.axis = 1.3)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/figure1_v2.pdf", width = 10, height = 5)

par(mfrow=c(1,2))
par(mar=c(4, 4, 1,0.1))
par(mgp=c(2.3,1,0))
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= expression('Fire size (km'^2*')'),ylim=c(0,200), cex.lab=1.4,cex.axis = 1.3)
text(0.3,195,"a)",pos=4, cex=1.4)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= "",ylim=c(0,200), cex.lab=1.4,cex.axis = 1.3)

text(0.3,195,"b)",pos=4, cex=1.4)
dev.off()

quartz_off_screen 
                3 

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_figure1_v2.pdf", width = 10, height = 5)

par(mfrow=c(1,2))
par(mar=c(4, 4, 1,0.1))
par(mgp=c(2.3,1,0))
boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= expression('Fire size (km'^2*')'),ylim=c(0,700), cex.lab=1.4,cex.axis = 1.3)
text(0.3,690,"a)",pos=4, cex=1.4)
boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("1","2","3","4","5"),xlab="Day since fire start",ylab= "",ylim=c(0,700), cex.lab=1.4,cex.axis = 1.3)

text(0.3,690,"b)",pos=4, cex=1.4)
dev.off() 

quartz_off_screen 
                3 

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/all_v3.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n',cex.axis=1.4, lty = 1, lwd = 2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
legend(x="topleft", legend=c("human","lightning"),col="black",lty = c(1,2),pch=1,bty = "n",lwd = 2,cex=1.5, pt.cex = 1)

dev.off() 

quartz_off_screen 
                3 

t.test(log(pro1$growth_km),log(pro1h$growth_km))

    Welch Two Sample t-test

data:  log(pro1$growth_km) and log(pro1h$growth_km)
t = 6.5083, df = 169.36, p-value = 8.265e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.492386 2.791866
sample estimates:
 mean of x  mean of y 
 1.3226886 -0.8194377 

t.test(log(pro2$growth_km),log(pro2h$growth_km))

    Welch Two Sample t-test

data:  log(pro2$growth_km) and log(pro2h$growth_km)
t = 6.9843, df = 134.4, p-value = 1.199e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 1.369607 2.451706
sample estimates:
mean of x mean of y 
2.7559564 0.8453003 

t.test(log(pro3$growth_km),log(pro3h$growth_km))

    Welch Two Sample t-test

data:  log(pro3$growth_km) and log(pro3h$growth_km)
t = 5.3694, df = 137.17, p-value = 3.286e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.9293526 2.0129234
sample estimates:
mean of x mean of y 
 3.186372  1.715234 

t.test(log(pro4$growth_km),log(pro4h$growth_km))

    Welch Two Sample t-test

data:  log(pro4$growth_km) and log(pro4h$growth_km)
t = 4.6217, df = 122.92, p-value = 9.467e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.6909157 1.7261008
sample estimates:
mean of x mean of y 
 3.513204  2.304696 

t.test(log(pro5$growth_km),log(pro5h$growth_km))

    Welch Two Sample t-test

data:  log(pro5$growth_km) and log(pro5h$growth_km)
t = 4.9328, df = 103.98, p-value = 3.089e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.9170094 2.1499736
sample estimates:
mean of x mean of y 
 3.895266  2.361774 

##################3 for western cordillera ecoregion ##################

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/north_v2.pdf", width = 4, height = 5)
par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2,cex.axis=1.4)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

null device 
          1 

for mediteranean california

pro1 = res[res$fire_day == 1 & res$human == 1 & (res$eco1 == 11),]
pro2 = res[res$fire_day == 2 & res$human == 1 & (res$eco1 == 11),]
pro3 = res[res$fire_day == 3 & res$human == 1 & (res$eco1 == 11),]
pro4 = res[res$fire_day == 4 & res$human == 1 & (res$eco1 == 11),]
pro5 = res[res$fire_day == 5 & res$human == 1 & (res$eco1 == 11),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & (res$eco1 == 11),]
pro2h = res[res$fire_day == 2 & res$human == 2 & (res$eco1 == 11),]
pro3h = res[res$fire_day == 3 & res$human == 2 & (res$eco1 == 11),]
pro4h = res[res$fire_day == 4 & res$human == 2 & (res$eco1 == 11),]
pro5h = res[res$fire_day == 5 & res$human == 2 & (res$eco1 == 11),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)


par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/med_v2.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="",cex.axis=1.4, xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                3 

for difference in autumn

pro1 = res[res$fire_day == 1 & res$human == 1 & res$month >8,]
pro2 = res[res$fire_day == 2 & res$human == 1 & res$month >8,]
pro3 = res[res$fire_day == 3 & res$human == 1 & res$month >8,]
pro4 = res[res$fire_day == 4 & res$human == 1 & res$month >8,]
pro5 = res[res$fire_day == 5 & res$human == 1 & res$month >8,]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & res$month >8,]
pro2h = res[res$fire_day == 2 & res$human == 2 & res$month >8,]
pro3h = res[res$fire_day == 3 & res$human == 2 & res$month >8,]
pro4h = res[res$fire_day == 4 & res$human == 2 & res$month >8,]
pro5h = res[res$fire_day == 5 & res$human == 2 & res$month >8,]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)


par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/autumn_v2.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,250),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n',cex.axis=1.4, lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                3 

for difference in summer

pro1 = res[res$fire_day == 1 & res$human == 1 & res$month <=8 ,]
pro2 = res[res$fire_day == 2 & res$human == 1 & res$month <=8,]
pro3 = res[res$fire_day == 3 & res$human == 1 & res$month <=8,]
pro4 = res[res$fire_day == 4 & res$human == 1 & res$month <=8,]
pro5 = res[res$fire_day == 5 & res$human == 1 & res$month <=8,]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & res$month <=8,]
pro2h = res[res$fire_day == 2 & res$human == 2 & res$month <=8,]
pro3h = res[res$fire_day == 3 & res$human == 2 & res$month <=8,]
pro4h = res[res$fire_day == 4 & res$human == 2 & res$month <=8,]
pro5h = res[res$fire_day == 5 & res$human == 2 & res$month <=8,]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)


par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/summer_v2.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n',cex.axis=1.4, lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                3 

for difference in summer in western cordillera

pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 6 | res$eco1 == 7),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

t.test(log(pro1$growth_km),log(pro1h$growth_km))

    Welch Two Sample t-test

data:  log(pro1$growth_km) and log(pro1h$growth_km)
t = 1.038, df = 53.703, p-value = 0.3039
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.4317907  1.3585976
sample estimates:
 mean of x  mean of y 
-0.2855046 -0.7489080 

t.test(log(pro2$growth_km),log(pro2h$growth_km))

    Welch Two Sample t-test

data:  log(pro2$growth_km) and log(pro2h$growth_km)
t = 3.6912, df = 82.447, p-value = 0.0003995
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.555548 1.854073
sample estimates:
mean of x mean of y 
 2.058908  0.854097 

t.test(log(pro3$growth_km),log(pro3h$growth_km))

    Welch Two Sample t-test

data:  log(pro3$growth_km) and log(pro3h$growth_km)
t = 2.2426, df = 64.873, p-value = 0.02835
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.07573427 1.30867147
sample estimates:
mean of x mean of y 
 2.535161  1.842958 

t.test(log(pro4$growth_km),log(pro4h$growth_km))

    Welch Two Sample t-test

data:  log(pro4$growth_km) and log(pro4h$growth_km)
t = 2.3845, df = 62.683, p-value = 0.02014
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.1141261 1.2961201
sample estimates:
mean of x mean of y 
 2.976315  2.271192 

t.test(log(pro5$growth_km),log(pro5h$growth_km))

    Welch Two Sample t-test

data:  log(pro5$growth_km) and log(pro5h$growth_km)
t = 2.3428, df = 40.216, p-value = 0.02417
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.1060498 1.4368888
sample estimates:
mean of x mean of y 
 3.327351  2.555881 

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)

par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/north_summer_v2.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="",cex.axis=1.4, xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                3 

for difference in autumn in western cordillera

pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro4h = res[res$fire_day == 4 & res$human == 2 & (res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 6 | res$eco1 == 7),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)

par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/north_autumn_v2.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,250),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n',cex.axis=1.4, lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                3 

for difference in summer in meditereanean

pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month <=8 & res$month >5 ) & (res$eco1 == 11),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

#t.test(log(pro1$growth_km),log(pro1h$growth_km))
#t.test(log(pro2$growth_km),log(pro2h$growth_km))
#t.test(log(pro3$growth_km),log(pro3h$growth_km))
#t.test(log(pro4$growth_km),log(pro4h$growth_km))
#t.test(log(pro5$growth_km),log(pro5h$growth_km))

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)

par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/med_summer_v2.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n',cex.axis=1.4, lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                3 

` ################## for difference in autumn in mediteranean ##################

pro1 = res[res$fire_day == 1 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro2 = res[res$fire_day == 2 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro3 = res[res$fire_day == 3 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro4 = res[res$fire_day == 4 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]
pro5 = res[res$fire_day == 5 & res$human == 1 & ( res$month >8 ) & (res$eco1 == 11),]

boxplot(pro1$growth_km,pro2$growth_km,pro3$growth_km,pro4$growth_km,pro5$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

pro1h = res[res$fire_day == 1 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro2h = res[res$fire_day == 2 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro3h = res[res$fire_day == 3 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro4h = res[res$fire_day == 4 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]
pro5h = res[res$fire_day == 5 & res$human == 2 & ( res$month >8 ) & (res$eco1 == 11),]

boxplot(pro1h$growth_km,pro2h$growth_km,pro3h$growth_km,pro4h$growth_km,pro5h$growth_km,names=c("day1","day2","day3","day4","day5"),xlab="",ylab="fire size (km2)",ylim=c(0,300), cex.lab=1.4,cex.axis = 1.3)

me[1] =mean(pro1$growth_km,na.omit=T)
me[2] =mean(pro2$growth_km,na.omit=T)
me[3] =mean(pro3$growth_km,na.omit=T)
me[4] =mean(pro4$growth_km,na.omit=T)
me[5] =mean(pro5$growth_km,na.omit=T)

me1[1] =mean(pro1h$growth_km,na.rm=T)
me1[2] =mean(pro2h$growth_km,na.rm=T)
me1[3] =mean(pro3h$growth_km,na.rm=T)
me1[4] =mean(pro4h$growth_km,na.rm=T)
me1[5] =mean(pro5h$growth_km,na.rm=T)

par(mgp=c(2.3,1,0))  

plot(me,type="o", ylim=c(0,150),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days)

points(me1,type="o",lty = 2, lwd = 2)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig4/med_autumn_v2.pdf", width = 4, height = 5)

par(mgp=c(2.3,1,0))  
plot(me,type="o", ylim=c(0,400),ylab=expression("fire size (km"^2*")"),xlab="", xaxt='n', lty = 1,cex.axis=1.4, lwd = 2,cex.lab=1.2)
axis(side=1, at=c(1:5),labels=days1,cex.axis=1.4)
points(me1,type="o",lty = 2, lwd = 2)
dev.off() 

quartz_off_screen 
                3 

how many days does it take to reach 75% burnt area

res75 = res[res$per_ba > 0.75,]
#peak_day = as.data.frame(aggregate(res75$fire_day, by = list(res75$firename,res75$cause), min))
#peak_day=subset(peak_day,x < 55)
#hi = hist(peak_day$x,prob =F, breaks= c(0:54), xlim=c(0,55), ylab="number of fires", xlab="days", cex.lab=1.4,cex.axis=1.3)

out1 = subset(res75,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 )


peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))
peak=as.data.frame(aggregate(res75$fire_day, by = list(res75$firename), min))
quantile(peak_day1$x,0.50,type=3) 

50% 
 10 

quantile(peak_day2$x,0.50,type=3) 

50% 
  3 

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,30))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/time_to_reach75.pdf",width=7,height=5)

fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,30))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  mediteranean  ###########
out1 = subset(res75,cause == 1 & res75$eco1 == 11)   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & res75$eco1 == 11)
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/med.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,15))

axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  north cal  ###########
out1 = subset(res75,cause == 1  & (res75$eco1 == 6 | res75$eco1 == 7))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14  & (res75$eco1 == 6 | res75$eco1 == 7))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/north.pdf",width=7,height=5)

fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,15))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  mediteranean  SUMMER ###########
out1 = subset(res75,cause == 1 & res75$eco1 == 11 & ( res75$month <=8 & res75$month >5 ))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & res75$eco1 == 11 & ( res75$month <=8 & res75$month >5 ))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/med_summer.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,10))

axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  north cal  summer ###########
out1 = subset(res75,cause == 1  & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month <=8 & res75$month >5 ))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14  & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month <=8 & res75$month >5 ))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/north_summer.pdf",width=7,height=5)

fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,12))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  mediteranean  autumn###########
out1 = subset(res75,cause == 1 & res75$eco1 == 11 & ( res75$month >8))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14 & res75$eco1 == 11 & ( res75$month >8 ))
#peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

#peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
#hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(0,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/med_autumn.pdf",width=7,height=5)
fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,5))

axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

########  north cal  autumn ###########
out1 = subset(res75,cause == 1  & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month >8))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res75,cause !=1 & cause != 14  & (res75$eco1 == 6 | res75$eco1 == 7)& ( res75$month >8 ))
peak_day1 = as.data.frame(aggregate(out1$fire_day, by = list(out1$firename), min))
peak_day2 = as.data.frame(aggregate(out2$fire_day, by = list(out2$firename), min))

peak_day1=subset(peak_day1,x < 56)
peak_day2=subset(peak_day2,x < 56)
hist.a =hist(peak_day1$x,breaks =c(0:55),plot=F)
hist.b =hist(peak_day2$x,breaks =c(0:55),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

pdf(file="/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/sup_fig5_time75/north_autumn.pdf",width=7,height=5)

fr = barplot(fg,xlab="Days after ignition",ylab="Number of fires",cex.lab=1.4,cex.axis = 1.3, xlim=c(1,65), ylim=c(0,5))
axis(1,c(0.7,5.5,11.5,17.5,23.5,29.5,35.5,41.5,47.5,53.5,59.5,65.5),labels=c(1,5,10,15,20,25,30,35,40,45,50,55),cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off()

quartz_off_screen 
                2 

res=res[res$ros_km>0,]
res_f = res[res$max_land == 1,]
res_p = res[res$max_land != 1,]

summary(lm(log(res$mean_frp)~log(res$ros_km)))

Call:
lm(formula = log(res$mean_frp) ~ log(res$ros_km))

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5266 -0.5420 -0.0497  0.4990  4.5702 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      1.74609    0.01869   93.41   <2e-16 ***
log(res$ros_km)  0.35186    0.01321   26.63   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8362 on 2249 degrees of freedom
  (547 observations deleted due to missingness)
Multiple R-squared:  0.2397,    Adjusted R-squared:  0.2394 
F-statistic: 709.1 on 1 and 2249 DF,  p-value: < 2.2e-16

#just show the plot here
plot(res_f$mean_frp~res_f$ros_km,log="xy",xlim=c(0.005,30),ylim=c(0.1,180),xaxt="n",ylab="mean FRP (MW)",xlab="Rate-of-Spread (km/day)", cex.lab=1.4,cex.axis = 1.3,col="darkgreen")

marks=c(0.01,0.1,1,10)
marks1=c(0.1,0.5,5,50)


pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/fig_FRP_ros_v3.pdf", width = 5, height = 5)

plot(res_f$mean_frp~res_f$ros_km,log="xy",xlim=c(0.005,30),ylim=c(0.1,180),xaxt="n",yaxt="n",ylab="mean FRP (MW)",xlab=expression('Rate-of-Spread (km d'^-1*')'), cex.lab=1.4,cex.axis = 1.3,col="darkgreen")
points(res_p$mean_frp~res_p$ros_km,col="orange")
axis(1,at=marks,labels=marks,cex.axis=1.4 )
axis(2,at=marks1,labels=marks1,cex.axis=1.4 )
legend( x="topleft",legend=c("Forest","Grass & shrub"),col=c("darkgreen","orange"),cex=1.2,pch=1,bty = "n")
dev.off()

quartz_off_screen 
                3 

FRP vs tree mortality

summary(lm(log(res$mean_frp)~res$mean_bas))

Call:
lm(formula = log(res$mean_frp) ~ res$mean_bas)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.3250 -0.5824  0.0428  0.5755  2.7943 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.853006   0.039957   21.35   <2e-16 ***
res$mean_bas 0.019447   0.001025   18.97   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8574 on 1727 degrees of freedom
  (1069 observations deleted due to missingness)
Multiple R-squared:  0.1725,    Adjusted R-squared:  0.172 
F-statistic: 359.9 on 1 and 1727 DF,  p-value: < 2.2e-16

summary(lm(data_forest$mean_BA_red~(log10(data_forest$ros_km)+I((log10(data_forest$ros_km))^2))))

Call:
lm(formula = data_forest$mean_BA_red ~ (log10(data_forest$ros_km) + 
    I((log10(data_forest$ros_km))^2)))

Residuals:
    Min      1Q  Median      3Q     Max 
-40.210 -11.114  -3.598   8.209  82.558 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       27.0251     0.5005   54.00   <2e-16 ***
log10(data_forest$ros_km)         25.2022     1.1296   22.31   <2e-16 ***
I((log10(data_forest$ros_km))^2)   8.3470     0.8597    9.71   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.51 on 1426 degrees of freedom
Multiple R-squared:  0.2908,    Adjusted R-squared:  0.2898 
F-statistic: 292.4 on 2 and 1426 DF,  p-value: < 2.2e-16

difference in fire size for first 5 days across california and both ecosystems

 
res$ros1 = res$max_ros+1


out1 = subset(res,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )

hum1 = out2[out2$fire_day ==1,]
hum2 = out2[out2$fire_day ==2,]
hum3 = out2[out2$fire_day ==3,]
hum4 = out2[out2$fire_day ==4,]
hum5 = out2[out2$fire_day ==5,]
lig1 = out1[out1$fire_day ==1,]
lig2 = out1[out1$fire_day ==2,]
lig3 = out1[out1$fire_day ==3,]
lig4 = out1[out1$fire_day ==4,]
lig5 = out1[out1$fire_day ==5,]
mean(hum1$growth)
mean(hum2$growth)
mean(hum3$growth)
mean(hum4$growth)
mean(hum5$growth)
mean(lig1$growth)
mean(lig2$growth)
mean(lig3$growth)
mean(lig4$growth)
mean(lig5$growth)
t.test(log10(hum1$growth),log10(lig1$growth))
t.test(log10(hum2$growth),log10(lig2$growth))
t.test(log10(hum3$growth),log10(lig3$growth))
t.test(log10(hum4$growth),log10(lig4$growth))
t.test(log10(hum5$growth),log10(lig5$growth))

hum1 = out2[out2$fire_day ==1 & out2$eco1==6,]
hum2 = out2[out2$fire_day ==2 & out2$eco1==6,]
hum3 = out2[out2$fire_day ==3 & out2$eco1==6,]
hum4 = out2[out2$fire_day ==4 & out2$eco1==6,]
hum5 = out2[out2$fire_day ==5 & out2$eco1==6,]
lig1 = out1[out1$fire_day ==1 & out1$eco1==6,]
lig2 = out1[out1$fire_day ==2 & out1$eco1==6,]
lig3 = out1[out1$fire_day ==3 & out1$eco1==6,]
lig4 = out1[out1$fire_day ==4 & out1$eco1==6,]
lig5 = out1[out1$fire_day ==5 & out1$eco1==6,]
mean(hum1$growth)
mean(hum2$growth)
mean(hum3$growth)
mean(hum4$growth)
mean(hum5$growth)
mean(lig1$growth, na.rm=T)
mean(lig2$growth, na.rm=T)
mean(lig3$growth, na.rm=T)
mean(lig4$growth, na.rm=T)
mean(lig5$growth, na.rm=T)
t.test(log10(hum1$growth),log10(lig1$growth))
t.test(log10(hum2$growth),log10(lig2$growth))
t.test(log10(hum3$growth),log10(lig3$growth))
t.test(log10(hum4$growth),log10(lig4$growth))
t.test(log10(hum5$growth),log10(lig5$growth))

hum1 = out2[out2$fire_day ==1 & out2$eco1==11,]
hum2 = out2[out2$fire_day ==2 & out2$eco1==11,]
hum3 = out2[out2$fire_day ==3 & out2$eco1==11,]
hum4 = out2[out2$fire_day ==4 & out2$eco1==11,]
hum5 = out2[out2$fire_day ==5 & out2$eco1==11,]
lig1 = out1[out1$fire_day ==1 & out1$eco1==11,]
lig2 = out1[out1$fire_day ==2 & out1$eco1==11,]
lig3 = out1[out1$fire_day ==3 & out1$eco1==11,]
lig4 = out1[out1$fire_day ==4 & out1$eco1==11,]
lig5 = out1[out1$fire_day ==5 & out1$eco1==11,]
mean(hum1$growth)
mean(hum2$growth)
mean(hum3$growth)
mean(hum4$growth)
mean(hum5$growth)
mean(lig1$growth, na.rm=T)
mean(lig2$growth, na.rm=T)
mean(lig3$growth, na.rm=T)
mean(lig4$growth, na.rm=T)
mean(lig5$growth, na.rm=T)
#t.test(log10(hum1$growth),log10(lig1$growth))
#t.test(log10(hum2$growth),log10(lig2$growth))
t.test(log10(hum3$growth),log10(lig3$growth))
t.test(log10(hum4$growth),log10(lig4$growth))
t.test(log10(hum5$growth),log10(lig5$growth))

out1 = subset(res,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )

mean(out1$ros_km,na.rm=T)

[1] 0.8405408

mean(out2$ros_km,na.rm=T)

[1] 1.82709

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7))

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

out1 = subset(res,cause == 1 & eco1 == 11)   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11)

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

############ autumn northern california
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7) & (month > 9))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7) & (month > 9))

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

############ autumn mediterean california
out1 = subset(res,cause == 1 & eco1 == 11 & (month > 9))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11 & (month > 9))

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

############ summer northern california
out1 = subset(res,cause == 1 & (eco1 == 6 | eco1 == 7) & (month > 5 & month<10))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & (eco1 == 6 | eco1 == 7) & (month > 5 & month<10))

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

############ summer mediterean california
out1 = subset(res,cause == 1 & eco1 == 11  & (month > 5 & month<10))   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 & eco1 == 11  & (month > 5 & month<10))

hist.a =hist(out1$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.01,0.05,0.1,0.25,0.5,1,2,5,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5,30.5,33.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

out1 = subset(res,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )


mean(out1$ros_km,na.rm=T)

[1] 0.8405408

mean(out2$ros_km,na.rm=T)

[1] 1.82709

#0,0.25,0.5,1,2,3,5,7,10,20,30
hist.a =hist(out1$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
hist.b =hist(out2$ros_km,breaks =c(0,0.5,1,2,3,5,7,10,20,30),plot=F)
fg = rbind(hist.a$counts,hist.b$counts)
dens = rbind((hist.a$counts/(sum(hist.a$counts)))*100,(hist.b$counts/(sum(hist.b$counts)))*100)

fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab="Number of fire days",cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

fr = barplot(dens, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab="% fire days",ylim=c(0,60),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/figure2_freq_v1.pdf", width = 6, height = 5)

fr = barplot(dens, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab="Fire days (%)",ylim=c(0,60),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off() 

quartz_off_screen 
                3 

###burnt area by rate-of-spread

data_test = data_s1[data_s1$max_land == 1,]
data_test1 = data_s1[data_s1$L1CODE == 6 |data_s1$L1CODE == 7 ,]
data_test2 = data_s1[data_s1$L1CODE == 11,]

data_test1 = data_s1[data_s1$human == 1 & (data_s1$L1CODE == 6 |data_s1$L1CODE == 7),]
data_test2 = data_s1[data_s1$human == 2 & (data_s1$L1CODE == 6 |data_s1$L1CODE == 7),]

data_test1 = data_s1[data_s1$human == 1 ,]
data_test2 = data_s1[data_s1$human == 2 ,]

fast = data_test1[data_test1$ros_km > 3.64,]
sum(fast$size,na.rm=T)/sum(data_test1$size,na.rm=T)

[1] 0.4992649

sum((data_test1$size*data_test1$ros_km),na.rm=T)/sum(data_test1$size,na.rm=T)

[1] 5.610449

fast = data_test2[data_test2$ros_km > 2.2,]
sum(fast$size,na.rm=T)/sum(data_test2$size,na.rm=T)

[1] 0.499989

sum((data_test2$size*data_test2$ros_km),na.rm=T)/sum(data_test2$size,na.rm=T)

[1] 3.532566

fast = data_s1[data_s1$ros_km > 1,]
fast_hum = fast[fast$human == 1,]

print("% BA and % of fire days fast fires > 1km/day")

[1] "% BA and % of fire days fast fires > 1km/day"

sum(fast$size)/sum(data_s1$size)

[1] 0.8359887

length(fast$size)/length(data_s1$size)

[1] 0.3371951

print("% BA  fast fires due to human ignition % of fire days human caused fast fires > 1km/day")

[1] "% BA  fast fires due to human ignition % of fire days human caused fast fires > 1km/day"

sum(fast_hum$size, na.rm=T)/sum(fast$size)

[1] 0.4634998

length(fast_hum$size)/length(fast$size)

[1] 0.4195298

all_min1 = data_s1[data_s1$nr_day != 1,] # remove first fire spread day from statistics

quan = quantile(data_s1$ros_km,0.9)
fast = data_s1[data_s1$ros_km > quan,]
slow = data_s1[data_s1$ros_km < quan,]
fast_hum = fast[fast$human == 1,]

print("fastest 10% fires cause xxx% of BA")

[1] "fastest 10% fires cause xxx% of BA"

sum(fast$size)/sum(all_min1$size)

[1] 0.5493523

length(fast$size)/length(all_min1$size)

[1] 0.1

print("mean tree mortality weighted by BA and just mean")

[1] "mean tree mortality weighted by BA and just mean"

sum((data_s1$mean_BA_red*data_s1$size))/(sum(data_s1$size))

[1] 49.03327

mean(data_s1$mean_BA_red)

[1] 25.64874

print("% BA due to human fires amoung fastest 10% fire days")

[1] "% BA due to human fires amoung fastest 10% fire days"

sum(fast_hum$size, na.rm=T)/sum(fast$size)

[1] 0.5143879

print("% fire number due to human fires amoung fastest 10% fire days")

[1] "% fire number due to human fires amoung fastest 10% fire days"

length(fast_hum$size)/length(fast$size)

[1] 0.5060976

print("% tree mortality fast fires weighthed and not")

[1] "% tree mortality fast fires weighthed and not"

sum((fast$mean_BA_red*fast$size))/(sum(fast$size))

[1] 60.02651

mean(fast$mean_BA_red)

[1] 52.79125

print("% tree mortality slow fires weighthed and not")

[1] "% tree mortality slow fires weighthed and not"

sum((slow$mean_BA_red*slow$size))/(sum(slow$size))

[1] 35.63221

mean(slow$mean_BA_red)

[1] 22.63291

print("tree mortality <0.5km and >2km")

[1] "tree mortality <0.5km and >2km"

fast1 = data_s1[data_s1$ros_km > 2,]
slow1 = data_s1[data_s1$ros_km < 0.5,]
mean(fast1$mean_BA_red,omit.na=T )

[1] 48.30526

mean(slow1$mean_BA_red,omit.na=T )

[1] 15.29343

# plot BA per rate of spread 
out1 = subset(data_s1,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(data_s1,cause !=1 & cause != 14 )

breaks =c(0,0.5,1,2,3,5,7,10,20,30)

tt=0
pp=0
for (i in 1:9){

  kr = out1[out1$ros_km >= breaks[i] & out1$ros_km < breaks[i+1],]
  kp = out2[out2$ros_km >= breaks[i] & out2$ros_km < breaks[i+1],]

tt[i]=sum(kr$size, na.rm=T)/1000000
pp[i]= sum(kp$size, na.rm=T)/1000000
}
sum(tt)

[1] 6781.34

sum(pp)

[1] 5612.935

fg = rbind(tt,pp)

par(mar=c(4, 5, 2,0.1))
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab=expression('burnt area (km'^2*')'),ylim=c(0,2000),cex.lab=1.4,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels= breaks ,cex.axis = 1.3)

legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")


tiff("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/BA_per_ros_hist.tif", width = 6, height = 5, units = 'in', res = 300)

par(mar=c(4, 5, 2,0.1))
fr = barplot(fg, beside=TRUE,xlab="Rate-of-Spread (km/day)",ylab=expression('burnt area (km'^2*')'),ylim=c(0,2000),cex.lab=1.4,cex.axis = 1.3)
#axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels= breaks ,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off() 

quartz_off_screen 
                3 

fg1=rbind((tt/sum(tt))*100,(pp/sum(pp))*100)

pdf("/Users/stijnhantson/Documents/Documents/articulos/en_proceso/VIIRS_ros/BA_per_ros_hist_percentage.pdf", width = 6, height = 5)
#par(mar=c(1, 1, 1,1))
fr = barplot(fg1, beside=TRUE,xlab=expression('Rate-of-Spread (km d'^-1*')'),ylab= "Burned area (%)",ylim=c(0,30),cex.lab=1.4,cex.axis = 1.3)

#axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels=hist.a$breaks,cex.axis = 1.3)
axis(1,at=c(0.5,3.5,6.5,9.5,12.5,15.5,18.5,21.5,24.5,27.5),labels= breaks ,cex.axis = 1.3)
legend("topright",legend = c("human","lightning"), fill=c("grey","black"),cex=1.4,bty = "n")
dev.off() 

quartz_off_screen 
                3 

20 fastest fires

res1 = res[res$ros_km > 10 & !is.na(res$ros_km),]

Error in res$ros_km : object of type 'closure' is not subsettable

are ROS the same for light & human under the same conditions

out1 = subset(res,cause == 1 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,cause !=1 & cause != 14 )

plot(out1$vpd,log(out1$ros_km))
points(out2$vpd,log(out2$ros_km),col="red")

summary(lm(out1$vpd~log(out1$ros_km+1),na.omit=T))
summary(lm(out2$vpd~log(out2$ros_km+1)))

analysis of the first day

load data

daily_res=read.table("/Users/stijnhantson/Documents/projects/VIIRS_ros/all_ignitions_V3.txt",header=T)

res=as.data.frame(daily_res)

res$bi =as.numeric(as.character(res$bi))
res$erc =as.numeric(as.character(res$erc))
res$etr =as.numeric(as.character(res$etr))
res$fm100 =as.numeric(as.character(res$fm100))
res$fm1000 =as.numeric(as.character(res$fm1000))
res$pet =as.numeric(as.character(res$pet))
res$pr =as.numeric(as.character(res$pr))
res$rmax =as.numeric(as.character(res$rmax))
res$rmin =as.numeric(as.character(res$rmin))
res$th =as.numeric(as.character(res$th))
res$tmmn =as.numeric(as.character(res$tmmn))
res$tmmx =as.numeric(as.character(res$tmmx))
res$vpd =as.numeric(as.character(res$vpd))
res$ws =as.numeric(as.character(res$ws))
res$vs =as.numeric(as.character(res$vs))
res$total_area =as.numeric(as.character(res$total_area))
res$max_land =as.numeric(as.character(res$max_land))
res$mean_land =as.numeric(as.character(res$mean_land))

res$biomass =as.numeric(as.character(res$biomass))

res = res[-1,]
res$human[res$cause ==1] =1
res$human[res$cause !=1 & res$cause !=14] =0

analysis

out1 = res[res$cause !=1 & res$cause != 14,] 
out2 = res[res$cause ==1,] 
length(out1$bi)

[1] 1196

length(out2$bi)

[1] 486

out1 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause !=1 & res$cause != 14)   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause ==1)   #1=lightning; 14=unknown; 7=arson
length(out1$bi)

[1] 335

length(out2$bi)

[1] 375

out1 = subset(res,eco1 == 11& res$cause !=1 & res$cause != 14)
out2 = subset(res,eco1 == 11 & res$cause ==1)
length(out1$bi)

[1] 784

length(out2$bi)

[1] 45

out1 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause !=1 & res$cause != 14 & res$mont > 5 & res$mont < 10 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause ==1& res$mont > 5 & res$mont < 10)   #1=lightning; 14=unknown; 7=arson
length(out1$bi)

[1] 248

length(out2$bi)

[1] 352

out1 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause !=1 & res$cause != 14 & res$mont < 6 & res$mont > 9 )   #1=lightning; 14=unknown; 7=arson
out2 = subset(res,(eco1 == 6 |eco1 == 7) & res$cause ==1& res$mont < 6 & res$mont > 9)   #1=lightning; 14=unknown; 7=arson

length(out1$bi)

[1] 0

length(out2$bi)

[1] 0

t.test(out1$bi,out2$bi)

Error in t.test.default(out1$bi, out2$bi) : not enough 'x' observations

analysis of MTBS temporal trend

library(raster)
library(rgeos)
mtbs_dir="/Users/stijnhantson/Documents/data/MTBS/DATA/"
med=shapefile("/Users/stijnhantson/Documents/data/veg_california/med_cal.shp")
north=shapefile("/Users/stijnhantson/Documents/data/veg_california/north_cal.shp")

year=1984
k=0
dnbr_all_med=0
rdnbr_all_med=0
dnbr_all_nor=0
rdnbr_all_nor=0

len_dnbr_med=0
len_rdnbr_med=0
len_dnbr_nor=0
len_rdnbr_nor=0
for (year in 1984:2017){
  print(year)
  k=k+1
mtbs_dir_year = paste(mtbs_dir,year,"/",sep="")
bndy_list = list.files(mtbs_dir_year, pattern = "burn_bndy.shp$", recursive = T, full.names=T)
shapefile_list <- lapply(bndy_list, shapefile)
fires <- do.call(rbind, shapefile_list)
fires=gUnaryUnion(fires)
fire_north = intersect(fires,north)
fire_med = intersect(fires,med)

dnbr = raster(paste(mtbs_dir,year,"_dnbr.tif",sep=""))
  dnbr[dnbr < -2000] <- NA
  rdnbr = dnbr = raster(paste(mtbs_dir,year,"_rdnbr.tif",sep=""))
  rdnbr[rdnbr < -2000] <- NA
  
dnbr_ext_med = extract(dnbr,fire_med)
dnbr_ext_nor = extract(dnbr,fire_north)

rdnbr_ext_med = extract(rdnbr,fire_med)
rdnbr_ext_nor = extract(rdnbr,fire_north)

dnbr_all_med[k]=mean(unlist(dnbr_ext_med),na.rm=T)
rdnbr_all_med[k]=mean(unlist(rdnbr_ext_med),na.rm=T)
dnbr_all_nor[k]=mean(unlist(dnbr_ext_nor),na.rm=T)
rdnbr_all_nor[k]=mean(unlist(rdnbr_ext_nor),na.rm=T)

len_dnbr_med[k] = length(unlist(dnbr_ext_med))
len_rdnbr_med[k] = length(unlist(rdnbr_ext_med))
len_dnbr_nor[k] = length(unlist(dnbr_ext_nor))
len_rdnbr_nor[k] = length(unlist(rdnbr_ext_nor))
removeTmpFiles(0)
gc()
}

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

```